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The  displacement  based  finite  element  model  of  a  general  third-order  beam  theory  is  developed  to  study 
the  quasi-static  behavior  of  viscoelastic  rectangular  orthotropic  beams.  The  mechanical  properties  are 
considered  to  be  linear  viscoelastic  in  nature  with  a  scope  to  undergo  von  Karman  nonlinear  geometric 
deformations.  A  differential  constitutive  law  is  developed  for  an  orthotropic  linear  viscoelastic  beam 
under  the  assumptions  of  plane-stress.  The  fully  discretized  finite  element  equations  are  obtained  by 
approximating  the  convolution  integrals  using  a  trapezoidal  rule.  A  two-point  recurrence  scheme  is 
developed  that  necessitates  storage  of  data  from  the  previous  time  step  only,  and  not  from  the  entire 
deformation  history.  Full  integration  is  used  to  evaluate  all  the  stiffness  terms  using  spectral /hp  lagrange 
polynomials.  The  Newton  iterative  scheme  is  employed  to  enhance  the  rate  of  convergence  of  the  non¬ 
linear  finite  element  equations.  Numerical  examples  are  presented  to  study  the  viscoelastic  phenomena 
like  creep,  cyclic  creep  and  recovery  for  thick  and  thin  beams  using  classical  mechanical  analogues  like 
generalized  n-parameter  Kelvin-Voigt  solids  and  Maxwell  solids. 
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1.  Introduction 

l.J.  Background 

In  the  practical  design  and  analysis  of  engineering  structures,  it 
becomes  very  important  to  predict  the  deflections,  strains  and 
stresses  to  prevent  a  catastrophic  failure.  If  these  are  made  of  vis¬ 
coelastic  materials,  then  it  becomes  critical  to  evaluate  the  re¬ 
sponse  of  the  structure  over  a  long  period.  Experimental  methods 
to  measure  these  are  often  costly,  time  consuming  and  sometimes 
not  possible.  The  theoretical  and  mathematical  background  of  vis¬ 
coelasticity  has  long  been  established  [1-7],  Most  of  the  analytical 
methods  use  the  correspondence  principle  [8],  to  analyze  visco¬ 
elastic  problems.  However,  this  method  is  restricted  to  very  limited 
problems  with  simple  geometry  and  loadings  for  which  explicit 
solutions  are  available. 

Hence,  for  practical  design  of  viscoelastic  structures  there  is  a 
need  for  numerical  methods  like  the  finite  element  method  [9], 
or  the  boundary  element  method  [10].  The  finite  element  method 
is  a  proven  technique  and  has  been  applied  to  static  and  dynamic 
problems  in  structural  mechanics.  Most  of  the  finite  element  ap¬ 
proaches  are  based  on  integral  transform  methods  [11],  or  super 
position  methods  [12].  However,  as  noted  by  Chen  and  Lin  [13], 
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these  methods  based  on  integral  transforms,  introduce  errors  into 
the  numerical  scheme,  due  to  the  approximate  nature  of  the  inver¬ 
sion  techniques.  They  presented  an  incremental  based  finite  ele¬ 
ment  technique  for  the  dynamic  response  of  viscoelastic  beams. 
It  avoids  the  integral  transformation,  and  thus  the  errors,  by 
replacing  the  creep  strain  increments  by  fictitious  body  forces. 

The  early  viscoelastic  finite  element  codes  are  based  on  history 
integral  forms,  as  in  Ref.  [14,15].  Finite  element  techniques  based 
on  these  methods  requires  the  storage  of  solution  vector  for  the  en¬ 
tire  deformation  history,  which  becomes  a  bottle  neck  when  the 
computational  memory  is  scarce.  Johnson  et  al.  [16]  developed  a 
technique  based  on  differential  constitutive  law;  Payette  and  Red¬ 
dy  [17]  and  Vallala  et  al.  [18]  developed  a  similar  technique  based 
on  a  two-point  recurrence  scheme.  In  the  finite  element  formula¬ 
tions  based  on  these,  the  data  storage  can  be  limited  to  the  last 
few  desired  history  steps.  These  formulations  make  use  of  mechan¬ 
ical  analogues  like  spring  and  dashpots  for  the  mathematical  mod¬ 
el,  to  predict  the  response  of  the  viscoelastic  structures.  The 
classical  generalized  models  of  Maxwell  solids,  Maxwell-Voigt  sol¬ 
ids  and  n-parameter  Kelvin-Voigt  solids  are  used  the  most. 

In  this  paper,  a  differential  constitutive  law  is  derived  based  on 
an  assumed  kinematic  field  (presented  in  the  next  section)  that 
necessitates  the  use  of  2D  plane-stress  constitutive  model.  The 
constitutive  relations  for  linear  anisotropic  viscoelastic  materials 
are  given  in  Ref.  [19,20].  We  make  use  of  fact  that  relaxation  mod¬ 
uli  for  a  linear  orthotropic  material  is  symmetric  as  verified  exper- 
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imentally  by  Halpin  and  Pagano  [20].  The  time  dependent  relaxa¬ 
tion  moduli  in  the  viscoelastic  constitutive  equation  is  expanded 
in  a  Prony  series  for  the  mechanical  analogue  models  used.  For 
more  on  the  constitutive  equations  and  the  number  of  independent 
constants  in  compliance/relaxation  tensor  see  Ref.  [21]. 

3.2.  Higher-order  beam  theories 


It  must  be  noted  that  the  assumed  form  of  displacement  field 
gives  rise  to  nonzero  strain  on  the  top  and  bottom  surfaces  of 
the  beam.  Also,  it  can  be  shown  that  all  other  third-order  and 
lower-order  theories  can  be  derived  from  this;  hence,  we  call  it  a 
general  third-order  beam  theory.  More  details  on  higher-order  shear 
deformable  theories  can  be  found  in  noted  works  by  Reddy  [22- 
24]. 


Beams  are  structural  elements  whose  length  is  large  compared 
to  their  cross-sectional  dimensions.  They  are  supported  at  few 
points  along  the  length,  and  subjected  to  forces  that  make  the 
structure  to  stretch  and  bend.  Theories  that  are  used  to  study  the 
response  of  beams  under  external  loads  are  obtained  by  reducing 
the  general  three-dimensional  elasticity  problem  through  a  series 
assumptions  concerning  the  kinematics  of  deformation  and  consti¬ 
tutive  behavior.  The  kinematic  assumptions  exploit  the  fact  that 
such  structures  do  not  experience  significant  strains  or  stresses 
associated  with  the  thickness  direction.  Thus,  the  solution  of  the 
three-dimensional  elasticity  problem  associated  with  a  beam 
structure  is  reformulated  in  terms  of  displacements  or  stresses, 
whose  form  is  presumed  on  the  basis  of  an  educated  guess  con¬ 
cerning  the  nature  of  the  deformation. 

Beam  theories  based  on  the  assumed  form  of  the  displacement 
field  are  most  popular.  In  these  theories,  the  displacements  are  ex¬ 
panded  in  increasing  powers  of  the  thickness  (or  height)  coordi¬ 
nate.  The  word  “order”  refers  to  the  power  of  the  thickness 
coordinate  in  the  power  series  expansion  of  the  displacement  field. 
The  cubic  expansion  of  the  displacement  field  is  optimal  because  it 
gives  quadratic  variation  of  transverse  shear  strain  and  stress,  and 
require  no  “shear  correction  factors”  compared  to  the  lower-order 
Timoshenko  beam  theory,  where  the  transverse  shear  strain  and 
stress  are  constant  through  the  beam  thickness. 

In  the  context  of  higher-order  theories,  there  is  no  beam  theory 
that  accounts  for  shear  deformation  while  not  requiring  shear  cor¬ 
rection  factors,  material  variation  through  the  beam  thickness,  and 
geometric  nonlinearity.  This  very  fact  motivated  the  present  study. 
The  objective  of  the  current  paper  is  to  develop  a  general  third-or¬ 
der  beam  theory  that  accounts  for  two-constituent  material  prop¬ 
erties  with  von  Karman  nonlinear  strains  to  capture  the  bending- 
extensional  coupling.  Hence  we  consider  displacement  field  of 
the  form 

iq  ( x ,  z,  t )  =  u(x,  t)  +  z9x(x ,  t)  +  z2<fx(x,  t)  +  z3il/x(x,  t ) 
u2{x,z ,  t)  =  0 

u3(x,z,t)  =  wo(x,t)+z0z(x,t)+z2<j>z(x,r)  (1.1) 


Then  the  nonzero  von  Karman  nonlinear  strains  are 
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2.  Constitutive  model 


As  shown  above,  the  assumed  displacement  field  makes  the 
strain  ezz  to  be  nonzero,  which  necessitates  us  to  consider  a  2D 
plane-stress  constitutive  model.  For  a  linear  orthotropic  material 
the  relation  between  second  Piola-Kirchhoff  stress  tensor,  denoted 
by  <t  and  reduced  strain  tensor  E  &  s  is  given  by 
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where  Q,j  are  the  time  independent,  plane-stress  reduced,  elastic 
coefficients  and  Qij  are  the  monotonically  decreasing  functions  of 
time  that  constitute  the  plane-stress  reduced  viscous  coefficients. 
The  elastic  plane-stress  coefficients  are  related  to  engineering 
material  constants  as 


1  *XZ  VZX  1  VXZ  VZX  1  *XZ  VZX 


q33(0)  =  ■ 


Ez(  0) 


Q55(0)  =  C*z(0) 


1  -  VxzVa  ’ 

the  time  dependent  viscous  stiffness  values  are 

4„<t 4„<t 


(2.2) 


Q33(t-s)  =  Ez(t  SJ  ,  Qss(t  -  s)  =  Cxz(t  -  s) 
1  —  vxzvzx 


(2.3) 


where  Ex(t )  and  Ez(t)  are  the  extensional  relaxation  moduli  and 
Gxz(t )  is  the  shear  relaxation  moduli.  Where  as  and  vzx  are  the 
major  and  minor  Poisson’s  ratios  of  the  beam.  The  specific  forms 
of  Ex(t),  Ez(t)  and  G^t)  will  depend  on  the  material  model  em¬ 
ployed.  In  this  study,  we  express  each  one  of  them  using  a  Prony 
series  of  order  n  as 


E(t)  =  E0  +  J2E,e'<,  G(t)  =  Go  +  ^C,eT 


(2.4) 


The  time  derivative  of  the  relaxation  moduli  can  be  expressed 
as 


y>a  —  sx+  +  z  (  2  4>x  + 


dth 

dx 
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ezz  =  6z  +  2  zif>z 

The  strain  field  can  be  expressed  as 


p  —  c(0)  70(1)  _[_  p(2)  ^3 p(3) 

C'XX  ~  c'xx  '  ^cxx  '  cxx  '  ^  cxx 

y*z  =  y{2+zy  {3)+z2y^ 

£zz  =  eg’  +  zeg1 


(1.2) 


(1.3) 


m  =  -±^Ee~i  G(t) 

1=1  zi 


(2.5) 


It  is  important  to  note  that  in  the  integral  constitutive  equations 
given  by  Eq.  (2.1)  we  assume  that  a  discontinuity  exists  in  the  re¬ 
sponse  only  at  f  =  0.  This  Prony  series  representation  of  the  visco¬ 
elastic  relaxation  moduli  is  critical  in  developing  the  recurrence 
scheme  and  to  implement  efficient  temporal  numerical  integration 
algorithms  of  the  viscoelastic  constitutive  equations. 


where 


3.  Weak  form  and  semi-discrete  element  model 


<o,  =  au  l/9w\2  m_d0x  (2)  =  d4  r(3)_d<l>x 

”  dx  2 \dxj  ’  **  dx’  ”  ax’  ”  dx 


dw 


de. 


(1.4) 


4^=2^  yS)=fls«+^,  yS>=2^+^,  (i.s) 


3.3.  Galerkin  weak  form 

To  derive  equations  of  motion  we  use  the  Hamilton’s  principle 
(see  Ref.  [23])  in  the  undeformed  configuration 
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[<5K~(<5l/  +  <5V)]  dt  =  0 


(3.1) 


where  SK  is  the  virtual  kinetic  energy,  SU  is  the  virtual  strain  energy, 
and  SV  is  the  virtual  work  done  by  external  forces.  Evaluating  each 
of  these  terms  and  performing  the  necessary  integration-by-parts 
with  respect  to  x  and  t,  we  obtain  the  following  weak  forms  of 
the  governing  equations  over  a  typical  finite  element  Qe  =  (xa,xb): 

0  =  J  ^mouSu  +  m^6xdu  +  m24>xdu+m3ij/xdu+Mi^^^-Fxdu^dx-  |  b 

(3.2) 

0  =  J  (m0wSw+mfBzdw+m2f>zSw+M^‘^^^+M^^^-Fz5w)dx 

-  [mS'^w+aOw]! 

(3.3) 


0= L  (mi  ii56x + m2^'5e* + m3  ^xSdx + m^,idex + dx 


(3.4) 


0  =  J  (^2ud$x+m36xd<l>x+m4<j)xdtl>ll  +  m5\j/xd<l>x+M$‘^‘l^+2M£<j<l>(SJdx 
-[Mg>^x]r 


(3.5) 


0=  I  (m3uSil/x+m48xS\l/:I+m5ij>xd\l/:l+m6ijix5tl/1[+M^^^^+3M^5il/x-fPs<l/^jdx 


(3.6) 


0  =  jf  b  (jn,w5Oz+m2dzdez  +  m3ij>zdOz+M(xl'l^+M(°>dOz-f^d0z'jdx 

~  [<^0Z 


(3.7) 


0  =  /  (m2'i',<5^z  +  +  mfaStfrz  +  MW  °^+  2M'Z'Z)S<I>Z dx 

(3.8) 


where  (Su,Sw,30x,Sc/>x  ,&\l/x,&dz,&4>z)  are  the  respective  virtual  varia¬ 
tions  that  can  be  viewed  as  the  test  functions,  and  m,  are  the  mass 
inertias  defined  by 

[m0,m,,m2,m3,m4,m5,m6]  =  /  p[1,z,z2,z3,z4,z5,z6]dA  (3.9) 

Ja 

The  various  internal  stress  resultants  are  defined  as 


a2Z,  aa\ z’dA 


(3.10) 


The  terms  Fx  and  Fz  due  to  the  external  loads  are  given  as 


Fx  =  /x(0) ,  Fz  =  /<°>  =>  /«  =  [ 1  zrfxdA.  f®  =  [Lz%dA  (3.11) 

Here  fx,fz  denote  the  distributed  axial  load  and  transverse  load 
about  the  y  axis  respectively.  The  variational  problem  for  the  gen¬ 
eral  third-order  beam  can  be  stated  as:  find  {u,w,Qx,4>xpj/x,  0Z,<^Z)  e 
H\G)  x  H\G)  x  H\Q)  x  H1  x  (£2)  H’(£2)  x  H1  (£2)  x  H\Q)  for  all 
( Su,dw,56x,5(l>x,dipx,Sez,S(l>z )  €  H1(£2)  x  H1  (£2)  x  H’(£2)  xH’x 


(£2)H1(£2)  x  H’(£2)  x  H’(£2)  such  that  the  equations  Eqs.  (3.2)- 
(3.8)  hold  true,  where  Hm(£ 2)  is  the  Sobolev  space  of  order  m  and 

12  =  [xa,Xk]. 


3.2.  Semi-discrete  finite  element  model 


Since  the  assumed  kinematic  displacement  requires  only  the 
continuity  of  the  primary  variables  across  the  element  boundaries 
and  not  its  derivatives,  i.e.,  C°-continuous,  we  use  the  following 
equal-order  interpolation  functions  for  all  variables 

[U(X,  t),W(X,  t),0x{X,  t),<f>x(X,  t),l//x{X,  t),ez{X,  t),  (fz(x,  t)] 
m 

=  l][uj(0,  Wj(t),  eXj(t),  (fxj(t),  xjjxj{t),  0g(t),  ^(WjM  (3.12) 
i=  i 


where  i/q  are  the  one-dimensional  nodal  spectral  interpolation 
functions.  The  nodal  expansion  in  the  interval  £2st=[-l,+l]  for  a 
typical  finite  element  is  defined  as 


m) 


(p)(p +  !)!„({,)({ -{j) 


(3.13) 


where  Lp  =  P°p°  are  the  Legendre  polynomials  of  order  p,  and  y  de¬ 
notes  the  location  of  the  roots  of  (£  -  1)({  +  l)Lp({)  =  0  in  the  inter¬ 
val  [-1.+1],  These  set  of  points  {§ Yff1  are  commonly  referred  to 
as  Gauss- Lobatto-Legendre  (GLL)  points.  Other  nodal  bases  like 
Chebyshev  second-kind  can  also  be  used.  The  location  of  nodes  in 
a  typical  master  element  coincides  with  the  roots  of  Legendre  poly¬ 
nomial  hence  the  basis  is  called  “spectral”  [25].  In  spectral  basis  the 
nodes  are  not  equally  spaced  within  the  canonical  interval  [— 1,+1], 
At  high  polynomial  orders  (five  and  above)  the  equi-spaced  interpo¬ 
lation  functions  exhibit  spurious  oscillations  at  the  ends  of  the 
interval,  this  is  called  Runge  effect.  This  impacts  the  accuracy  and 
reliability  of  finite  element  formulation.  The  spectral  nodal  basis 
does  not  suffer  from  this  and  also  the  Kronecker  delta  property, 
i.e.,  i j/j  ({,)  =  Sjj  is  not  compromised. 

In  Fig.  la  we  can  see  Runge  oscillations  for  equally  spaced  La¬ 
grange  basis  near  the  ends  of  the  interval  and  in  Fig.  lb  the  cardi¬ 
nality  condition  of  Kronecker  delta  is  obvious.  In  generating  the 
above  plot  and  in  our  code  an  analytically  less  complex  and  com¬ 
putationally  more  stable  form  of  Eq.  (3.13)  is  used  to  generate 
the  nodal  basis.  The  equation  is  shown  below 


P+1 

n 

j=i 

jVi 


Uc-£,-) 

te  -  (j) 


(3.14) 


Multi-dimensional  spectral  interpolation  functions  can  be  ob¬ 
tained  from  the  tensor  product  of  the  above  one-dimensional 
equation  [9,26].  Substituting  Eq.  (3.12)  into  the  weak  forms  in 
Eqs.  (3.2)-(3.8),  yields  the  semi-discrete  finite  element  model  of 
the  beam  element.  The  quasi-static  finite  element  equations  can 
be  expressed  as  (and  are  given  explicitly  in  Appendix  A) 


ds  =  F 


(3.15) 


4.  Fully  discretized  finite  element  equations 

4.1.  Time  discretization  using  recurrence  formulation 

In  order  to  derive  fully  discretized  finite  element  equations, 
we  start  with  the  partitioning  of  the  time  interval  [0.  T]  c  R 
(region  of  interest)  into  set  of  N  non-overlapping  subintervals 
such  that 


3762 


V.  Vallala  et  al. / Composite  Structures  94  (2012)  3759-3768 


Lagrange  Interpolation  Functions 


Fig.  1.  Graphs  of  (a)  equi-spaced  and  (b)  spectral  lagrange  interpolation  functions  for  polynomial  order  of  p  =  11. 


N 


[0,  T|  =  U 

k=  1 


4+i  ] 


(4.1) 


The  final  solution  is  obtained  by  repeatedly  solving  an  initial  va¬ 
lue  problem  within  each  subregion  [4,tfc+1],  with  the  known  values 
of  solution  at  f  =  tk  as  initial  conditions.  From  Eq.  (3.15)  it  is  clear 
that  the  semi-discrete  finite  element  equations  have  contributions 
from  elastic  and  viscous  parts  of  the  constitutive  relations.  The 
elastic  part  is  simple  and  straight  forward,  however,  the  contribu¬ 
tion  from  viscous  part  is  in  the  form  of  convolution  integrals,  hence 
the  full  discretization  of  these  is  not  a  trivial  task.  In  order  to  solve 
the  problem  in  each  subinterval,  we  can  approximate  these  convo¬ 
lution  integrals  using  two-point  (trapezoidal  rule)  or  three-point 
(Simpson’s  rule)  formulas.  But  a  direct  temporal  integration  from 
here,  results  in  a  computationally  unattractive  solution  procedure 
which  requires  the  storage  of  the  entire  deformation  history.  It  be¬ 
comes  a  bottle  neck  for  storing  these  when  the  computational 
memory  is  scarce.  Also,  when  the  total  number  of  time  steps  N  is 
large,  much  of  the  computational  time  expended  to  get  the  solu¬ 
tion  at  a  subinterval,  goes  into  the  evaluation  of  the  convolution 
integrals. 

To  circumvent  these  issues,  we  develop  a  recurrence  scheme  for 
two-point  (trapezoidal  rule)  formula  that  can  be  used  to  approxi¬ 
mate  the  convolution  integrals  with  in  each  subinterval.  The 
two-point  recurrence  scheme  requires  the  storage  of  the  general¬ 
ized  displacements  and  a  set  of  internal  variables  evaluated  at 
the  Gauss  points,  from  the  previous  time  step  only.  A  similar 
three-point  Simpson's  scheme,  which  requires  the  storage  from 
last  two  time  steps,  is  used  in  Ref.  [27].  Using  these  ideas,  the  con¬ 
volution  integral  appearing  in  Eq.  (3.15)  can  be  expressed  as 

/*t n  _  1  rh+ i 

/  KA(s)  ds  =  J2  KA(s)  ds  (4.2) 

Jo  fc=l  Jtk 

In  order  to  develop  the  recurrence  formulation  the  following 
multiplicative  decomposition  of  the  relaxation  moduli  [28],  is  used. 
These  equations  hold  true  as  the  relaxation  moduli  can  be  ex¬ 
pressed  in  terms  of  Prony  series  within  each  subinterval. 

£(4+i  -  s)  =  ^VAtk/Ti i£|(tk  -  s) 

/=  1 

G(4+i  -  s)  =  Q(t*  -  s)  (4.3) 

/=  1 


where  A  tk  =  tk+i  -  t k.  Using  the  above,  Eq.  (4.2)  can  be  expressed  in 
index  notion  at  an  arbitrary  time  step  t  =  ts  as 

s- 1  ptk+x  ___  n  NGP 

X,(ts)  =  /  KijAjis)  ds  S  J]5>mX!m(ts)  (4.4) 

k=  1  1=1  m= 1 


where  Einstein’s  summation  convention  on  repeated  indices  is  im¬ 
plied.  As  noted  previously,  Gauss  quadrature  is  employed  in  evalu¬ 
ation  of  Ky,  resulting  in  the  summation  over  m  (where  NGP  is  the 
number  of  Gauss  points).  The  quantity  Xjm  assumes  the  following 
possible  forms  for  extensional  and  shear  moduli 


Xjm(4)  =  e  lfx;m(ts_,)“ 


At,  i  E, 


e  (4-5) 


ms  l  A  f-  c 

xj'"(fs)  =  e^rx;m(ts_1)-%i5 

z  L, 


*'  r(4-,)+/im(4)  (4.6) 


Note  to  get  the  above  expressions,  we  replaced  convolution 
integrals  with  in  each  subinterval  with  a  two-point  trapezoidal 
rule.  Also,  the  specific  forms  of  xm  and  J)m(G)  depend  on  compo¬ 
nents  of  I<jj.  In  Eq.  (4.4),  even  though  there  are  (s  -  1)  time  steps 
of  k  to  get  to  time  t=  ts,  the  above  equations  just  need  the  values 
of  solution  {A(ts)}  and  internal  variables  X[m(ts_,)  from  k  =  ts  and 
k  =  ts_i .  There  is  no  need  to  store  these  values  for  all  the  (s  -  1) 
time  steps.  Thus  the  above  equations  represent  recurrence  formu¬ 
las  in  terms  of  the  internal  variables  Xjm(ts)  with  X[m(t1  =  0)  =  0. 
Using  Eqs.  (4.5),  (4.6)  results  in  the  following  quasi-static  fully-dis¬ 
cretized  equations  for  generalized  displacements  at  the  current 
time  step  (the  components  are  presented  in  Appendix  B) 

[X]s{A(ts)}  =  {F}s-{Q}s  (4.7) 


5.  Numerical  examples 

In  this  paper,  we  make  use  of  higher-order  spectral  interpola¬ 
tion  functions  without  resorting  to  any  selective  or  reduced  inte¬ 
gration  techniques.  The  nonlinear  finite  element  equations  are 
linearized  using  Newton-Raphson’s  procedure  (see  Appendix  B) 
and  the  equations  are  solved  using  an  iterative  scheme.  Since  a 
nonlinear  beam  becomes  stiff  with  load,  the  total  load  is  divided 
into  smaller  load  steps,  with  the  solution  of  each  step  being  used 
for  the  next.  For  all  the  problems  we  use  five  load  steps  with  a 
maximum  of  10  nonlinear  iterations  at  each  step.  All  problems  in 
this  study  converged  within  five  iterations  for  a  tolerance  of 
e=  10~6.  The  computational  domain  is  reduced  by  taking  advan¬ 
tage  of  the  symmetry  of  the  beam  about  x  =  1/ 2.  The  boundary  con¬ 
ditions  considered  in  the  analysis  are. 

(1)  Hinged  at  both  ends: 

w(  0)  =  0,  u(L/2)  =  0,  0X(L/ 2)  =  0,  4>x{L/2) 

=  0,  WL/2)=0  (5.1) 

(2)  Pinned  at  both  ends: 

u(  0)  =  0,  w(0)  =  0.  u(L/2)  =  0,  ex(L/2) 

=  0,  <j>x(L/2)  =  0,  UL/2)  =  0 


(5.2) 
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Table  1 

Viscoelastic  moduli. 


Eo 

E, 

E2 

e3 

e4 

Es 

Ee 

Et 

Es 

E9 
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s  5 


53 

<D 
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2  - 


205.7818  Ksi 
43.1773  Ksi 
9.2291  Ksi 
22.9546  Ksi 
26.2647  Ksi 
34.6298  Ksi 
40.3221  Ksi 
47.5275  Ksi 
46.8108  Ksi 
58.6945  Ksi 


(1) 


(1) 


9.1955  x  10-’  s 
9.8120  x  10°  s 
9.5268  x  101  s 
9.4318  x  102  s 
9.2066  x  103s 
8.9974  x  104s 
8.6852  x  105  s 
8.5142  x  106  s 
7.7396  x  107  s 
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Fig.  3.  Quasi-static  maximum  vertical  deflection  wmax,  of  hinged-hinged  beam 
under  harmonic  time-dependent  transverse  loading  q(t). 


Bakker  [29],  for  a  glassy  amorphous  polymer  material  known  as 
PMMA.  The  Prony  series  parameters  for  the  viscoelastic  relaxation 
modulus  given  in  Table  1  were  calculated  by  Payette  and  Reddy 
[17],  from  the  published  compliance  parameters  in  Ref.  [29]. 
Although  the  finite  element  formulation  places  no  restriction  on 
the  relationship  between  £(t)  and  G(t),  for  the  present  analysis 
we  adopt  the  approach  taken  by  Chen  [30],  and  assume  that  the 
shear  and  relaxation  moduli  are  related  as  G(t)  =  E(t)/ 2(1 +v), 
where  v  is  Poisson’s  ratio  of  the  material,  which  is  assumed  to  be 


Time  1 ,  [s] 

Fig.  2.  Quasi-static  maximum  vertical  deflection  wmax,  of  viscoelastic  beam  under 
uniform  distributed  load  q. 

Table  2 

Quasi-static  finite  element  results  for  the  maximum  deflection  wmax  of  a  viscoelastic 
beam  under  uniform  distributed  load  q  with  three  different  sets  of  boundary 
conditions. 


time-independent  and  equal  to  v  =  0.4  [31]. 

A  viscoelastic  beam  of  uniform  cross  section  lxl  in.,  and 
length  L  =  100  in.,  with  material  properties  given  in  Table  1  is  used 
for  the  analysis.  At  t  =  0  the  beam  is  subjected  to  a  time  invariant 
uniform  vertical  distributed  load  q  =  0.25  lbf/in.  A  constant  time 
step  At  =  1.0  s  is  employed  with  a  total  simulation  time  of 
1800  s.  Graphical  results  are  presented  in  Fig.  2  for  a  beam  discret¬ 
ized  with  two  finite  elements  with  5th  order  spectral  interpolation 


Time,  t  Fliigge 
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8.5429 

400 

8.6827 
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8.9886 

1600 

9.0271 
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9.0612 

At  =  0.1  At  =  0.5 


7.2995 

7.2995 

8.5457 

8.5648 

8.6856 

8.7052 

8.7710 

8.7910 

8.8394 

8.8597 

8.8976 

8.9182 

8.9478 

8.9687 

8.9917 

9.0127 

9.0302 

9.0514 

9.0644 

9.0857 

At  =  1 .0  At  =  2.0 


7.2995 

7.2995 

8.6237 

8.8512 

8.7661 

9.0014 

8.8531 

9.0931 

8.9228 

9.1666 

8.9820 

9.2291 

9.0333 

9.2832 

9.0780 

9.3304 

9.1172 

9.3719 

9.1520 

9.4086 

(3)  Clamped  at  both  ends: 

u(  0)  =  0.  w(0)  =  0,  0X(O)  =  0,  t/>x(  0)  =  0,  i/q(0)  =  0 

u(L/2)  =  0.  0*(L/2)  =  0,  cj>x(L/2)  =  0,  i/c,(L/2)  =  0 

(5.3) 


5.1.  Quasi-static  analysis  of  thin  beams 

In  the  first  example,  we  compare  the  results  for  a  thin  isotropic 
beam  ( L/h  >  20)  in  x  and  z  directions  with  the  results  presented  in 
Ref.  [18],  for  a  n-parameter  Kelvin-Voigt  solid.  The  viscoelastic 
material  model  is  based  on  the  experimental  findings  of  Lai  and 
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Fig.  4.  Quasi-static  maximum  vertical  deflection  wmax,  of  hinged-hinged  elastic  and 
viscoelastic  beams  under  time-dependent  transverse  loading  q{t). 
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Fig.  5.  Quasi-static  maximum  vertical  deflection  wmax,  of  hinged-hinged  beam 
under  heaviside  time-dependent  transverse  loading  q(t). 


functions.  As  expected,  the  deflection  steadily  increases,  and  then 
approaches  a  zero  slope  of  equilibrium  or  a  long-time  constant  va¬ 
lue.  This  behavior  is  called  creep  under  constant  load.  Also  at  t  =  0, 
the  results  coincide  with  the  instantaneous  elastic  solution,  where 
Young’s  modulus  is  given  as  E  =  535.4 Ksi. 

For  the  hinged-hinged  beam  configuration,  the  vertical  deflec¬ 
tion  is  compared  with  the  following  exact  solution  by  Fliigge  [8], 
for  a  geometrically  linear  viscoelastic  beam  based  on  the  Timo¬ 
shenko  beam  theory 


w0(L/2,t) 


5q0L4  [  8(1  +  v)  fh\ 
384 D2  5  Ks  \LJ 


D(t) 


(5.4) 


where  D(f)  is  the  creep  compliance  and  I(s  is  the  shear  correction 
factor.  The  results  are  given  in  Table  2  with  different  values  of  time 
steps  At. 

In  the  second  example,  we  study  the  cyclic  creep  response  of  a 
thin  orthotropic  beam  modeled  as  a  Maxwell  Solid.  The  values  of 
elastic  moduli  taken  from  Johnson  et  al.  [16],  are  £x  =  38.6  GPa, 
Ez  =  8.27  GPa,  Cxz  =  4.14  GPa,  =  0.26,  vzx  =  v^EJEx  and  the  Prony 
series  is  taken  as 

P(t)  =  1.0  +  0.01755e-°ooolt  +  0.000257e-°olt 

+  0.072014e-°3162t  (5.5) 

The  time  dependent  viscous  relaxation  moduli  is  taken  as 
Q[t)  =  QP(t).  A  beam  with  dimensions,  length  L  =  0.6  m,  height 
h  =  0.02  m,  and  base  b  =  0.01  m  is  used  for  analysis.  The  beam  is 
simply  supported  and  is  subjected  to  the  following  harmonic  dis¬ 
tributed  load  on  the  surface. 


5.2.  Quasi-static  analysis  of  thick  beams 

In  this  example,  we  present  numerical  solutions  for  creep  re¬ 
sponse  of  thick  (L/h  <  20)  orthotropic  beam  modeled  as  a  Maxwell 
solid.  The  elastic  moduli  and  the  Prony  series  are  the  same  as  in  the 
above  example.  The  dimensions  of  the  beam  are  length  L  =  0.1  m, 
height  h  =  0.02  m,  and  base  b  =  0.01  m.  The  beam  is  simply  sup¬ 
ported  and  is  subjected  to  a  vertical  distributed  load  on  the  top  sur¬ 
face.  The  load  is  ramped  from  a  value  of  0.0-1. 0  MPa  in  0.05  s  and 
is  maintained  constant  for  the  rest  of  the  time.  A  constant  time  step 
At  =  0.01  s  is  employed  with  a  total  simulation  time  of  5.0  s.  The 
beam  is  discretized  with  1 0  finite  elements  with  4th  order  spectral 
interpolation  functions.  The  creep  response  is  shown  in  Fig.  4,  as 
expected  both  the  elastic  and  viscous  response  converge  over  time. 

Next,  we  investigate  the  elastic  creep  recovery  behavior  of  the 
constitutive  model  using  the  following  time  dependent  load  for 
the  same  beam  as  above.  An  important  characteristic  is  that  the 
beam  should  eventually  return  to  its  original  configuration  once 
the  loads  are  removed.  To  demonstrate  that  the  present  finite  ele¬ 
ment  model  captures  this  effect,  we  consider  a  hinged-hinged 
beam  subjected  to  the  below  quasi-static  transverse  load  of  inten¬ 
sity  q0 

q(t)  =  q0  jw(t)  ~  T{j}l_  [(£ -  «T)H(t  -  ar)  -  (t  ~  P? )H(t  -  /It)]  j 

(5.7) 

where  H(t)  is  the  Heaviside  function,  and  we  take  q0  =  1.0  MPa  and 
x  =  1800  s.  The  parameters  0  ^  a  ft  ^  1  are  constants.  The  Eq.  (5.7) 
represents  a  load  function  that  is  constant  in  0  <  t  <  ax  and  then  lin¬ 
early  decreases  to  zero  from  t  =  act  to  f  =  /It.  For  f  >  fix,  the  load  is 
maintained  at  zero  (see  the  inset  in  Fig.  5).  A  constant  time  step 
At  =  0.01  s  is  employed  with  a  total  simulation  time  of  18  s  and 
the  beam  is  discretized  with  10  finite  elements  with  4th  order 
spectral  interpolation  functions.  In  Fig.  5  we  present  the  numerical 
results  for  various  values  of  a  and  p.  As  expected,  each  one  of  the 
curves  in  the  figure  follow  a  path  of  delayed  recovery  from  f  =  at 
to  t  =  p x  and  to  their  original  configurations  as  time  t  gradually 
tends  to  infinity  with  the  applied  loads  being  removed. 

6.  Conclusions 

In  this  paper  a  fully-discretized  finite  element  model  for  an 
orthotropic,  linear  viscoelastic  beam  based  on  a  general  third-order 
beam  theory  has  been  developed.  The  2-D  plane-stress  constitutive 
model  is  used  for  the  viscoelastic  formulation.  The  beam  consid¬ 
ered  is  capable  of  undergoing  moderate  rotations  and  small  strains 
in  the  sense  of  the  von  Karman  geometric  nonlinear  strains.  The  as¬ 
sumed  displacement  field  allows  for  a  better  bending-extensional 
coupling  and  the  use  of  C°-continuous  functions  for  all  the  primary 
variables,  thus  simplifying  the  implementation.  A  two-point  recur¬ 
rence  scheme  is  developed  such  that  history  from  the  last  time  step 
needs  to  be  stored.  Various  numerical  examples  have  been  in¬ 
cluded  to  demonstrate  the  capabilities  of  the  developed  finite  ele¬ 
ment  model. 


q(t)  =  0,  t  0 

<?(£■)  =  <Jo  [sin  +</>)],  t>0  (5.6) 

where  q0  =  0.01  MPa,  x  =  0.04  and  </>  =  0.  The  beam  is  discretized 
using  40  finite  elements  with  3rd  order  spectral  interpolation  func¬ 
tions.  A  constant  time  step  At  =  0.001  s  is  employed  with  a  total 
simulation  time  of  1  s.  In  Fig.  3  we  present  the  values  of  maximum 
vertical  deflection  of  the  beam.  The  oscillations  in  the  maximum 
and  minimum  values  demonstrate  the  cyclic  creep  behavior  of  the 
beam. 


Acknowledgments 

The  research  reported  herein  was  carried  out  while  the  first  and 
third  authors  were  supported  by  the  MURI09  Project  from  the  Air 
Force  Office  of  Scientific  Research  under  Grant  FA9550-09-1-0686. 

Appendix  A.  Semi-discretized  equations 

In  this  Appendix  we  present  explicitly  the  components  of  the 
semi-discretized  finite  element  model. 
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Kf^pD^t-sf^pdx,  Kf  =  £D„Q13(f-5)f^x 
Kfpy^-^dx,  Kr  =  £‘D2Q„(t-s)gf§dx 

Kf=*pD2Q55(t-s)3^jdx 

Kf  =  £D„Q13(t-s)Jf^dx,Kf%£iD2Q55(t-S)f§dx 

Kf  =  J [  DoQnV-sW^dx 

Kf  =  pDAu(t-s}^d^dx+D0(ls5(t-s)Mjdx, 

Kf  =  p  D4an(t-s)d^a^dX+D2(l5S(t-s)3Mjdx 
Kf  =  p  2D2Qn(t-s)d^<kjdx+D2(k5(t-sWi8^dx 
Kf„pD2(ln(t-s)^dX 

^P^-sP^^dx 

Kf=  p  D4Qn(t-s)(ppdx+2D2(l55(t-s)2^dx 

Kf  =  pD2Q_13(t-s)^jdx+2D2Clss(t-s)^dx 

Kf  =  jf  3D2Q55(t-s)^dx 

Kf  =  p  D4Qn(t-s)^  pdx  +  3D2  Q55(t- s)^dx, 

Kf  =  p  D6Q_u  (t- s)  dp  pdx+ 9D4Q55  (t  -  s^iijdx 
—■  rXb  ■  dil/-  •  dil/- 

Kf  =  J  2D4Q,3(t-s)^dx+3D4a55(t-s)^-^dx 
kf  =  pD2Q55(t-s)2^jdx+D2y3(t-s)^i,pdx 

Kf  =  p  D2  Q55  ( t  ■ -  s)  ^  ^  dx + D0Q33  (t  ■ -  S)  Mjdx 

*H°^-^dx 

kf  =  pD2Q55(t-s)tp^dx+2D2Q_,3(t-s)^Cpdx 
~  /*xb  .  dil/-  •  dil/- 

Kf  =  J  D4Q55(t-s)3^dx+2D4Q,3(t-s)^-^dx 

Kf  =  p  D4Q55(t-  s)  ^  pdx + 4D2  Q33  (t  -  s)<l/i>l/jdx  (A.5 


)  ^  -gfdx+D0Q33  (t  -  fOi/^dx 
)^dx 


All  other  stiffness  coefficients  are  zero. 


Appendix  B.  Fully-discretized  equations  and  tangent  matrix 

B.J.  Fully-discretized  finite  element  equations 

The  additional  matrices  introduced  in  the  fully-discretized  form 
of  the  finite  element  equations  can  be  expressed  as 
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D.(a„(0)  +  %lQ„(0))f 

0l(Qlim+^i4lim)jgaa„ 


‘  Do  (q,3  (0)  +  ^  Qi3  (0))  ^  *jdx 

Vq„<o>+^«„(o,)~* 


•D„(Q,1(0)  +  %lQ„(0))l(|:)2f  ^*  +  Da(a,,(0)  +  %lQ5!(0))|i§dx 

6Do^55(0)  +^=iQ55(0))  ^dx, 

b  n  ( rt  ,n,  AtN_i  A  ,n  A  d^j  dlj/j 

D2(Qll(0)+  —  an(0)j-^^dx, 

6  D2  ^0.55  (0)  +  ^1  Qss(O))  3  ^  fjdx 

tDo(Q.^O)+^Q.a(P))^*idx, 


'D2(a55(o)+^a55(o))f^dx 

1  Do  (q55(0)  +  ^-055(0))  <ki  ^dx, 


6  D2  (q„(0)  +  ^Qu(O))  ^  ^dx  +  D0  (q55(0)  +  ^ziQ55(0))  ^dx, 

‘  D4  (q_u( 0)  +  ^=1  Qn (0))  d^9^dx  +  D2  £55(0)  +^^Q55(0)^)  3 M}dx 
‘  2  D2  (q13(0)  +  ^(2,3(0))  ^  fjdx  +  D2  (Q55(  0)  +  ^Qs5(0))  d£dx 


6D2^11(0)+^lQ11(0))^^dx, 
6D4fQn(0)  +  A^lQ„(0)N)^ig/dx  +  2D2 


^  ^dx  +  2D2  (q55(0)  +  ^QssWW^dx, 


1  =  ^D2(q13(0)  +^lQ13(0)^lAJ.dx  +  2D2(^Q55(0)  +^lQ55(0)^i^'dx 

=  £  3 D2  ^Q.55(0)  +  ^-Q55(0))  *, ^dx, 

=  £  D4  £,  (0)  +  ^Qn(O))  ||^dx  +  3 D2  (q55(0)  +  ^|z1q55(0))  ^ dx , 
=  £d6£u(0)  +^Qn(0))^  ^dx  +  9D4(Q55(0)  +^|ziQ55(0))^.dx, 
=  £  2D4£3(0)  +^<2,3(0))  ^jdx  +  3D4 (^Q55(0)  +  A£^Q55(0)£i^dx 


1  Do  (q13  (0)  +  Qi3(0  )£i^dx, 

‘o.(a,3(0)+^Q„(0))’f*,|!* 
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K?.5  : 
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K77 
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:  £  D2  ^Q.55(0)  +  ^(255(0)  j  ^'Vjdx  +  2D2  (q,3(0)  +  ^Q13(0))^^  dx 
=  £  D4  ^Q.55(0)  +  ^Q55(0))3  ^%dx  +  2D4  (q13(0)  +  ^=1q,3(0))  ^dx 
=  I ^  D4  (W 0)  +  ^QssfO))  ^  ^ dx  +  4D2 ( Q33( 0)  +^=1033(0)N) M }dx 


(B-2) 


All  other  terms  are  zeros.  The  specific  forms  of  the  viscoelastic 
force  vector  Q  in  Eq.  (4.7)  and  the  history  terms  JX|m(tN)  appearing 
in  it  are  not  provided  due  to  space  constraints;  the  derivation  of 
these  components  is  similar  to  the  one  shown  in  Ref.  [18]. 

B.2.  Tangent  stiffness  matrix 


The  fully  discretized  finite  element  equations  are  nonlinear  due 
to  inclusion  of  the  von  Karman  strains  in  the  formulation.  For  our 
formulation  we  solve  the  equations  iteratively  using  Newton- 
Raphson  linearization  procedure  [32],  The  linearized  equations 
are  of  the  form 


{A(r)}s  =  {A(r-’>}s 

-  m;1  (r-1)],{Atr-,)}I  +  {Fir-%  -  {Q(r-1,}s)  (B.3) 

where  {A(r)}s  represents  the  solution  at  the  fth  iteration  and  time 
t=  ts.  The  tangent  stiffness  matrix  [T]s  in  Eq.  (B.3)  is  defined  using 
Einstein’s  summation  notation  as 


Tij  =  J2K<J  + 

1 


dKim 

dAj 


Am  + 


dQj 

dAj 


(B.4) 


All  quantities  in  Eq.  (B.4)  comprising  the  tangent  stiffness  ma¬ 
trix  are  formulated  using  the  solution  from  (r  -  1  )'th  iteration,  it 
is  important  to  note  that  all  the  partial  derivatives  are  taken  with 
respect  to  the  solution  of  the  current  time  step.  Applying  the  New¬ 
ton’s  method  to  the  fully-discretized  beam  equations  results  in  the 
following  components  of  tangent  stiffness  matrix 
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